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Abstract 

Recently developed quantum algorithms suggest that in principle, quantum computers can solve problems 
such as simulation of physical systems more efficiently than classical computers. Much remains to be done 
to implement these conceptual ideas into actual quantum computers. As a small-scale demonstration of 
their capability, we simulate a simple many-fermion problem, the Fano-Anderson model, using liquid state 
Nuclear Magnetic Resonance (NMR). We carefully designed our experiment so that the resource require- 
ment would scale up polynomially with the size of the quantum system to be simulated. The experimental 
results allow us to assess the limits of the degree of quantum control attained in these kinds of experiments. 
The simulation of other physical systems, with different particle statistics, is also discussed. 

PACS numbers: 03.67.-a, 05.30.-d, 76.60.-k, 03.65.Yz 



* Corresponding author. Email: camille@iqc.ca 



1 



I. INTRODUCTION 
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Quantum mechanical systems provide new resources to solve problems which are difficult to 
solve on classical computers. If we had a large quantum computer today, we could break crypto- 
lic codes [I], perform a variety of search algorithms flQ], estimate eigenvalues of operators 
], or simulate quantum systems [6]. In particular, the latter would enable a better understand- 
ing of the quantum world by enabling analyses of complex chemical reactions or demonstrating 
new states of matter. However, questions like What are the physical quantum states that can be 
reached efficiently? or What kind of physical processes can be efficiently simulated on a quantum 
computer? still remain open. 

Since Richard P. Feynman conjectured that an arbitrary discrete quantum system may be sim- 
ulated by any other |aj, the simulation of quantum phenomena became a fundamental problem 
that a quantum computer, i.e., a universally controlled quantum system, may potentially solve in a 
more efficient way than a classical computer. The basic idea is to imitate the evolution of a phys- 
ical system by cleverly controlling the evolution of the quantum computer. Quantum simulation 
is the process of faithfully imitating a physical phenomenon using a quantum computer. Although 
Feynman's illuminating conjecture seems appealing, it was only recently proved generally valid 
BZllMISlllOllllI]- Experimentally demostrating that one has universal control and thus can quantum 
imitate an arbitrary physical process constitutes an extremely challenging enterprise. 

It is important to notice that the efficiencies of quantum simulating the evolution of a physical 
system and of obtaining the sought-after information about a physical property must be established 
sep arately in most cases. A demonstration that evolution can be simulated efficiently kJ 
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1211 . that is, can be simulated with polynomial resources as a function of problem size, is 



in general insufficient for showing that the desired property (e.g., the ground state energy of a 
given Hamiltonian) can be obtained efficiently also. In general, the exponentially large Hilbert 
space that characterizes those physical systems and the inherent quantum parallelism of a quantum 
computer are insufficient for showing that an algorithm for quantum computation efficiently solves 
a problem. We pointed out in JsMll]] that in a quantum computation, it is necessary to demonstrate 
that in addition to maintaining adequate accuracy (noise, approximations, and statistical error 
control) one also has to demonstrate the polynomial scaling of the three main steps of a simulation, 
initialization, propagation and measurement. 

Some quantum processes can be simulated very well and efficiently on classical computers. 
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Simulating quantum phenomena using stochastic approaches reduces the problem to quadratures, 
which are multidimensional integrals that can be computed using Monte Carlo techniques. In gen- 
eral, the complexity of deterministic ^-dimensional integration is of order £- N / a (i.e., exponential 
in N), where e < 1 is some stipulated error and a quantifies the smoothness of the integrand. On 
the other hand, the expected complexity of Monte Carlo integration is of order e~ 2 , and hence 
independent of N and a (assuming that the variance of the integrand is finite). The reason for 
introducing these statistical techniques was to overcome the exponential complexity of determin- 
istic approaches such as the Lanczos method Jl^l . Realistic models of liquid or solid 4 He have 
been simulated to experimentally measured precision for a few years ifl^L Recently developed 
loop-cluster algorithms allow highly efficient and informative simulation of many quantum spin 
models of magnetism Il5 1. 

An important class of problems for which classical computers have major difficulties is the 
simulation of interacting fermionic systems (almost all large-scale simulations of fermions are 
done by the Monte Carlo method). In fact, as noted in JslQ], Feynman and others prior to him 
intuited this difficulty. Unless an approximation is made, the various quantum Monte Carlo algo- 
rithms must inevitably sample from a multivariate distribution P that has regions of phase space 
where it is negative that are comparable to regions where it is positive (because the state func- 
tion belongs to the totally antisymmetric representation of the permutation group). In general, the 
nodal hyper-surface P = separating the regions is unknown (an exception being when symme- 
try considerations alone determine it), making it impossible to solve the problem by independently 
sampling from each region where P has a definite sign. The sign problem is prohibitive on a classi- 
cal computer because it results in the variance of measured quantities growing exponentially with 
the number of degrees of freedom of the system. Still other applications require sampling from a 
complex-valued distribution P. This occurs, for example, if the simulation is done as a function 
of real Minkowski time or if time-reversal symmetry is broken. In previous work JsKk]], we have 
discussed how certain sign problems can be overcome using quantum network algorithms. 

In this paper we describe how quantum simulation of many-body problems can be realized in 
liquid state NMR Quantum Information Processors (QIPs) [16]. The constituents of the system 
may represent particles with arbitrary exchange statistics and generalized Pauli exclusion principle 
(such as fermions obeying Fermi statistics), spins, etc. In particular, we show how to efficiently im- 
itate a resonant impurity (localized state) scattering process in a metal (which is made of fermions), 
using the nuclear spins of a trans-crotonic acid molecule. This problem is physically modeled by 
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a Fano-Anderson Hamiltonian ||8|]. Our results demonstrate that the universal control achieved 
by the liquid state NMR QIPs enables efficient simulation of some fermionic (and other particle 
statistics) systems, providing relevant information about the particular phenomenon or system of 
study 111711 . In particular, we show how the spectrum of the Fano-Anderson Hamiltonian can be 
determined. 

The paper is organized in the following way: In Sec. |n|we introduce the conventional model 
of quantum computation and use it to describe the physics of the liquid state NMR setting as 
a universal quantum simulator. In Sec. [Hi] we show quantum algorithms for obtaining relevant 
physical properties of quantum systems satisfying different particle statistics, by mapping their 
algebras of operators into the spin- 1/2 algebra (conventional model). In Sec. [IV] we introduce the 
fermionic Fano-Anderson model, and show how to simulate it in the liquid state NMR device. Its 
experimental implementation as well as the results, and the conclusions are described in Sec. |yj 
and Sec.|yjJ respectively. 



H. QUANTUM INFORMATION PROCESSING WITH LIQUID STATE NMR METHODS 



In this section we introduce liquid state NMR quantum information processing methods, em- 
phasizing the fact that they can be mathematically described in terms of Pauli (spin- 1/2) operators 
1 18|]. A more detailed description of such methods can be found in I i 

In the conventional model of quantum computation the fundamental unit of information is the 
quantum bit or qubit. A qubit's pure state, |a) = a|0) + 6|1) (with a,b £ C and \a\ 2 + \b\ 2 = 1), 
is a linear superposition of the logical states |0) and |1), and can be represented by the state of a 
two-level quantum system such as a spin- 1/2. Similarly, a pure state of a register of N qubits is 
represented as \ip) = ^2 n=0 a n \ n ), where \n) is a product of states of each qubit in the logical 
basis, e.g., its binary representation (|0) = |00 ■ ■ • 0), |1) = |00 • • ■ 01), |2) = |00 - • • 10), etc.), and 
Yln=Q l \ a n\ 2 = 1 ( a n £ C). A quantum register can also be in a probabilistic mixture of pure 
states, i.e., a mixed state, which is described by a density matrix p = J2 s P s P s > w ^ m P s = \^l J s)('4 ) s\ 
representing the state of the register in the pure state |^ s ), with probability p s . Every density 
operator can be written as a sum of products of the Pauli spin- 1/2 operators a 3 a (a = x, y, z, and 
j — [1, • • • , N]) and the identity operators P acting on the j-th qubit of the register 111 611 . 

The Pauli operators can also be used to describe any unitary operation acting on the state of the 
register. In particular, every unitary operation can be decomposed in terms of single-qubit rotations 
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mum, defining a universal set of elementary gates. In Fig. [l]we show the quantum circuit 
representation of these basic operations. 

Finally, in the conventional model of quantum computation the measurement is assumed to be 
projective and is described by projectors that can be expanded in terms of Pauli operators. 

Liquid-state NMR methods allow us to physically implement a slightly different version of the 
conventional model of quantum computation, with respect to the initial state and the measurement 
process. In this set-up the quantum register is represented by the average state of the nuclear spin- 
1/2 of an ensemble of identical molecules. Since all molecules are equivalent, in the following 
analysis we will first consider only one of them. The spin state of each nucleus (qubit) of a single 
molecule is manipulated using resonant radio-frequency magnetic pulses (RF pulses). 

The molecule is placed in a strong magnetic field B(z) ~ 10 Tesla, so that the spin of the j-th 
nucleus precesses at its (Larmor) frequency Vj (Fig. In the frame rotating with the j-th spin, 
its qubit state can then be rotated by sending RF pulses in the x-y plane at the resonant frequency 
v r ~ Vj. If the duration of this pulse is At, the corresponding evolution operator in the rotating 
frame is lllfll 



where A is the amplitude of the RF-pulse and is its phase in the x-y plane (H = 1). Then one 
can induce single spin rotations I21I along any axis in the x-y plane by adjusting At and tp. 

Single-qubit rotations around the z-axis can be implemented with no experimental imperfection 
or physical duration simply by changing the phase of the abstract rotating frame we are working 
with. We have then to keep track of all these phase changes with respect to a reference phase 
associated with the spectrometer. Nevertheless, these phase tracking calculations are linear with 
respect to the number of pulses and spins, and can be efficiently done on a classical computer. 
Together with the rotations along the x- or y-axis, the z -rotations can generate any single qubit 
rotation on the Bloch sphere. 

On the other hand, the spin-spin interactions present in the molecule allow us to perform two- 
qubit gates and achieve universal control. To first order in perturbation, this interaction (called the 
J-coupling), has the form 



—iA(cos(ip)a J x +s'm(ip)ay )At 



(1) 



Hj,k = —a J z a z , 



(2) 
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where j, k denote the corresponding pair of qubits and Jjk is their coupling strength. Under typ- 
ical NMR operating conditions, these interaction terms are small enough to be neglected when 
performing single-qubit rotations with RF pulses of short duration . Nevertheless, between two 
pulses they are driving the evolution of the system. By cleverly designing a pulse sequence, i.e., 
a succession of pulses and free evolution periods, one can easily apply two-qubit gates on the 
state of the system. Indeed, the so-called refocusing techniques' principle consists of performing 
an arbitrary Ising gate by flipping one of the coupled spins (n -pulse), as shown in Fig. |3] The 
interaction evolutions before and after the refocusing pulse compensate leading to the effective 
evolution 



where the effective coupling strength a = Jjk(Ati — At 2 ) is being determined by the difference 
between the durations At\ and At 2 - 

We have so far described a quantum register as consisting of nuclei of a single molecule. How- 
ever, liquid state NMR uses an ensemble of about 10 23 molecules in a solution maintained at room 
temperature (~ 300A'). For typical values of the magnetic field, this thermal state is extremely 
mixed. Clearly, this is not the usual state in which we initialize a quantum computation since qubits 
are nearly randomly mixed. Nevertheless, known NMR methods iflfl can be used to prepare the 
so-called pseudo-pure state (p pp ) 
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Ppp = nN I + ePpure, (4) 



where p pure is a density operator that describes a pure state and e is a small real constant (i.e., e 
decays exponentially with N). 

Under the action of any unitary transformation U this state evolves as 

plf = Up pp U^ = ^r-I + Ue Ppure U j . (5) 

The first term in Eq. |5]did not change because the identity operator is invariant under any unitary 
transformation. Therefore, performing quantum computation on the ensemble is equivalent to 
performing quantum computation over the initial state represented only by p pure . 

After the quantum computation is performed, we measure the orthogonal components of the 
sample polarization in the x-y plane, M x = Tr(p f p 'p al J2f=i a l)> an ^ M y = Tr(p f '" al J2f=i a% y )- Note 
that the invariant component of p f p 'p al does not contribute to the signal since Tr(Ia 3 x ) = 0. Since 
the polarization of each single spin, M 3 X = Tr(p f p 'p al cr^) and M 3 y = Tr(p f p 'p al a^), precesses at its own 



Larmor frequency Uj, a Fourier transformation of the temporal recording (called FID, for Free 
Induction Decay) of the total magnetization needs to be performed. By doing so, we obtain the 
expectation value of the polarization of each spin (averaged over all molecules in the sample). 

Summarizing, a liquid state NMR setting allows us to initialize a register of qubits in a pseudo- 
pure state, apply any unitary transformation to this state by sending controlled RF pulses or by 
free interaction periods, and measure the expectation value of some quantum observables (i.e., the 
spin polarization). Hence, these systems can be used as quantum information processors (QIPs). 



III. SIMULATION OF PHYSICAL SYSTEMS 

Richard R Feynman f6] described a quantum computer as a universal reversible device gov- 
erned by the laws of quantum physics and capable of exactly simulating any physical system. 
Although he analyzed the problem of simulating physics assuming that every finite quantum me- 
chanical system can be imitated exactly by another one (e.g., a set of qubits) |7J], he was unsure 
whether this statement remained valid for the simulation of fermionic systems. 

In this section we describe how to obtain information about physical properties of any quantum 
many -body system (fermionic, bosonic, anyonic, etc.) by using a set of qubits (spin- 1/2) controlled 
by NMR techniques. A more complete description of these methods based on the existence of 
one-to-one mappings between the algebras used to describe the system to be simulated and the 
quantum computer Iq , [ill I24I . as well as indirect measurement algorithms [8], can be found in 
previous works JsL flO . Esfl . 

In this work we are interested in the measurement of correlation functions of the form 

G(t) = (4>\U(t)\4>), (6) 

where U(t) is any time (or other continuous parameter) dependent unitary operator, using indirect 
measurement techniques |8]. In addition to the qubits used to represent the physical system to be 
simulated (i.e., the system of qubits), an extra qubit called ancilla is required (Fig. EJ). This qubit 
will be used as a probe to scan the properties of the system of qubits. It has to be initialized in 
the superposition state |+) a = ^4=^ by applying the Hadamard gate ] to the polarized state 
|0) a . Then, it interacts with the system of qubits, initially in the state \(p), through a controlled 
unitary operation U' 1 ^ = |0) a (0| <g> I + |l) a (l| ® U{t). After this interaction, we can show 181] that 
G(t) = (2cx+) = {a 3 x + i<Jy); that means we get the desired result by measuring the expectation 



7 



values of the ancilla qubit observables a*, and a*. 

Using the same techniques we can determine the spectrum of an observable Q when choosing 
U(t) = e" 1 ® 1 . Figure |5] depicts this algorithm licll . Since the initial state can always be written as 
a linear combination of eigenstates of Q, that is, |0) = ^TnlVVi)' w i m l^n) m e eigenstates of Q 

n 

having eigenvalues A n , and 7 n complex coefficients, a measurement on the polarization of the an- 
cilla qubit gives (2cr+(t)) = |7„| 2 e _;iAn *. Having the time-dependent function S(t) = (2a\(t)) 

n 

for a discrete set of values U, the eigenvalues A n can in principle be obtained by performing a dis- 
crete Fourier transform (DFT) lip! . Note that the determination of each single value S(ti) requires 
a different experiment. 

The eigenvalues A n denote the spectrum of a system Hamiltonian H when replacing Q — > H. 
In this case, the operation U' < a can be efficiently implemented However, methods 

for finding an initial state with an overlap 7 n that does not vanish exponentially with increasing 
system size, are in general not known. This issue arises, for example, when trying to obtain the 
spectrum of the two-dimensional Hubbard model approaching the thermodynamic limit I loll^ l. 

Nevertheless, the same basic procedure can be used when interested in obtaining dynamical 
correlation functions of the form G(t) = (4>\T^AiTBj\<f)) (i.e., U(t) = T^A-iTBj in Eq. 0, 
where T = e~ %m is the time evolution operator of a time-independent Hamiltonian H, and Aj, Bj 
are unitary operators. In Fig. |6]we show the circuit for an algorithm capable of obtaining these 
correlation functions after some simplifications lip! . The evolution has three different steps: First, 
we perform a controlled operation B' 1 ^ = |0) a (0| <8> / + |l) a (l| ® Bj. Second, we perform the T 
operation on the system, and third, a controlled operation A^ a = |0) a (0| <g> A\ + |l) a (l| <S> I. Spatial 
correlation functions can also be obtained when replacing the operator T by the space translation 
operator. Again, this algorithm can be performed efficiently whenever the initial state \<p) can be 
prepared efficiently. 

The algorithm described above can be easily implemented with liquid- state NMR methods, 
since the result of the simulation is encoded in the expectation values of single qubit observables. 
So far, the algorithm applies only to the simulation of systems described in terms of Pauli oper- 
ators, such as spin-1/2 systems. However, other systems with different particle statistics can also 
be simulated with these algorithms after mapping their operator algebras onto the Pauli spin-1/2 
algebra lEKlill^l]. In the next section we introduce the Fano- Anderson model, a simple fermionic 
system, and show how to simulate it on a liquid-state NMR QIP using these methods. 
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IV. THE FANO-ANDERSON MODEL 



The quantum simulation of the one-dimensional fermionic Fano-Anderson model provides a 
starting point for simulations of quantum systems with different kinds of particle statistics. 

The one-dimensional fermionic Fano-Anderson model consists of an n-sites ring with an im- 
purity in the center (see Fig. HJ), where spinless fermions can hop between nearest-neighbors sites 
with hopping matrix element (overlap integral) r, or between a site and the impurity with matrix 
element V/ -Jn. Taking the single-particle energy of a fermion in the impurity to be e, and consid- 
ering the translational invariance of the system, the Fano-Anderson Hamiltonian can be written in 
the wave vector representation as [§] 

n-1 

H = E £ *A c h + ^ b + V (cl b + tfc k0 ), (7) 
1=0 

where the fermionic operators c£ and fet (c fe; and b) create (destroy) a spinless fermion in the 
conduction mode k\ and in the impurity, respectively. Here, the wave vectors are ki = — (I = 
[0, .., n — 1]) and the energies per mode are = — 2r cos ki. 

In this form, the Hamiltonian in Eq. |7]is almost diagonal and can be exactly solved: There are 
no interactions between electrons in different modes ki, except for the mode k , which interacts 
with the impurity. Therefore, the relevant physics comes from this latter interaction, and its spec- 
trum can be exactly obtained by diagonalizing a 2 x 2 Hermitian matrix, regardless of n and the 
number of fermions in the ring N e . Nevertheless, its simulation in a liquid-state NMR QIP is the 
first step in quantum simulations of quantum many -body problems. 

In order to use the algorithms presented in Sec. |nij and to successfully simulate this system 
in an NMR QIP, we first need to map the fermionic operators onto the spin- 1/2 (Pauli) operators. 



This is done by use of the following Jordan-Wigner transformation 112411 



b = a l _ tf = al 



c k 



(8) 



In this language, a logical state |0j) (with |0) = |f) in the usual spin-1/2 notation) corresponds 
to having a spinless fermion in either the impurity, if j = 1, or in the mode kj-2, otherwise. The 
fermionic vacuum state |vac) (i.e., the state with no fermions) maps onto |vac) = 1 1 ! 1 2 ■ • ■ 
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(= I X 1 4-2 • • • |n+i))- As an example, Fig.|7]shows the mapping of a particular fermionic state for 
n — 4. 

Some dynamical properties of this model can be obtained using the quantum algorithms de- 
scribed in Sec. UTTJ Here, we are primarily interested in obtaining the probability amplitude of 
having a fermion in mode k at time t, if initially (t = 0) the quantum state is the Fermi sea state 

N e -l 

with N e fermions; that is, |FS) = [ c' k |vac). This probability is given by the modulus square 
of the following dynamical correlation function: 

G(t) = (FS|6(t)& t (0)|FS), (9) 

where b{t) = T t 6(0)T, T = e~ im is the time evolution operator, and ^(O) = tf. Basically, G(t) 
is the overlap between the quantum state &t(0)|FS), which does not evolve, and the state &t(i)|FS), 
which does not vanish unless the evolved state T|FS) already contains a fermion in the impurity 
site ((fet(t)) 2 = (V(0)) 2 = 0). In terms of spin-1/2 operators (see Eq. [HJ>, this correlation function 
reduces to a two-qubit problem J^: 

G(t) = (0|TVT^|0) , (10) 
where T = e~ tHt is an evolution operator arising from the interaction terms in Eq. |7J with 

tt ^ 1 I S ko 2 | V („ 1 J2, | n n 

2 ' ~Y z * x + v y> ' ^ ' 

and |0) = | I1O2) in the logical basis (i.e., the initial state with one fermion in the k mode). 

In order to use the quantum circuit shown in Fig.|6l all operators in Eq. [I0|must be unitary. Us- 
ing the symmetries of H, such as the global ir/2 z-rotation that maps (a x , a J y ) — > (o y , —cr J x ), leav- 
ing the state \4>) invariant (up to a phase factor), we obtain (0|TVj.Tcri|0) = (cf)\T^ayTa x \(f)) = 
and {(j)\T^a x Ta x \(f)) = {(f>\T^GyTay\4>). Then, Eq. [TO]can be written in terms of unitary operators 
as 

G(t) = {<j>\e iat ale- iSt a x \<j>). (12) 

Figure |9] shows the quantum circuit used to obtain G(t). It is derived from Fig.|6]by making the 
following identifications: T — > e~ lHt , — > o\, and Bj — > a\. As we can see, the corresponding 
controlled operations A' ^ and B' 1 ^ transform into the well-known controlled-not (CNOT) gates. 
All the unitary operations appearing in Fig. were decomposed into elementary NMR gates 
(single qubit rotations and Ising interactions). In particular, the decomposition of e~ tHt can be 
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found in Ref. H8H. We obtain 

e~ im = Ue-'^e-'^U* , (13) 

where A 1(2) = \{E T VA 2 + V 2 ), with E = ^f*&, and A = The unitary operator U is 

decomposed as (Fig. |9]> 

U = e^ve^^e^&^e^v^^ (14) 

with the parameter 9 satisfying cos 6=1/ \A + S 2 , and 5 = ( A + VA 2 + V 2 ) /V. 

The CNOT gates Al°> a and B^ can also be decomposed into elementary gates, obtaining 
Al°>» = |0) a (0| ®crl + |l) a (l| <8> I = e^e^^e-^e-^^e-^ and Bl 1 ^ = |0) a (0| ® 
/ + 1 1 ) a ( 1 1 <g> o\ = e l % ax e~ t % a *^e~ t %°ve t * <r * a *e~ t % <T * (up to a phase factor). In this way, we can 
proceed to simulate the circuit of Fig. |9] and obtain G(i) in an NMR QIP by applying the appro- 
priate RF pulses (Sec. |n|). Only three qubits are required for its simulation (Fig. |9j): The ancilla 
qubit a, one qubit representing the impurity site (qubit-1), and one qubit representing the k mode 
(qubit-2). 

We are also interested in obtaining the spectrum of the Hamiltonian H of Eq. |7J For this 
purpose we used the algorithm shown in Fig. replacing Q — > H. In particular, when n = 1 (one 
site plus the impurity), Eq. |7] reduces to H = e+ ^ k ° + H, with H defined in Eq. [TT]in terms of 
Pauli operators. In this case, the two eigenvalues Aj (i = 1, 2) of the one-particle subspace can be 
extracted from the correlation function (SeclTTlh 

S(t) = (<f)\e- im \<j)} = e- i(e+£fe o)*(0| e -^|0), (15) 

which is equal to the polarization of the ancilla qubit after the algorithm of Fig. [5] is performed. 
Since |0) = |li0 2 ) = IJ.1T2) is not an eigenstate of H, it has a non-zero overlap with the two 
one-particle eigenstates, called |lPj) (see Appendix lAl). 

Again, the operator e lHa '^ t ^ 2 (Fig. |3]) needs to be decomposed into elementary gates for its 
implementation in an NMR QIP. Noticing that [cr|, H] = [cr|, U] = 0, we obtain 

e iHa%t/2 _ jj e i\ 1 a 1 z alt/2 e iX 2 a 2 z alt/2jj J i e i(e+e kQ )alt/2 ^ 

where the unitary operator U is decomposed as in Eq. [JU Figure [TO] shows the corresponding 
circuit in terms of elementary gates. Again, qubits 1 and 2 represent the impurity site and the k 
mode, respectively, a denotes the ancilla qubit. Since the idea is to perform a DFT on the results 
obtained from the measurement (see Appendix[A]), we need to apply this circuit for several values 

oft (Sec. unt- 
il 



V. EXPERIMENTAL IMPLEMENTATION 



A. Experimental protocol 

For the experimental simulation of the fermionic Fano-Anderson model, we used an NMR QIP 
based on a solution of trans-crotonic acid and methanol dissolved in acetone. This setting has 
been described in Ref. [27]. Once the state of the 3 equivalent protons in the methyl group of 
the trans-crotonic acid molecule is projected onto the spin- 1/2 subspace I27I . this molecule can be 
used as a seven-qubit register (see Fig. ITTT) . Methanol is used to perform RF-power selection and 
accurately calibrate the RF pulses. 

Two important characteristics of a molecule used for an NMR QIP are: (i) the accuracy of the 
control and (ii) the number of elementary gates we can perform within the relevant decoherence 
time of the system. The accuracy of control in trans-crotonic acid has been determined in Ref. 
1 2811 . using an error-correcting code as a benchmark. The current experiment can be considered as 
another exploration of the accuracy of control, in this case examining how well we can implement 
the necessary evolutions when simulating quantum systems with NMR techniques. 

In liquid-state NMR the main source of decoherence is the relaxation of the transversal polar- 
ization of the sample due to the loss of coherence between molecules. In our setting, the relevant 
times of this process, called T%, are in the range from several hundreds of milliseconds to more 
than one second, for the different nuclei. These times fix the maximum number of elementary gates 
that can be applied to the quantum register without lossing coherence. Indeed, a lower bound of 
the pulse duration to induce a rotation on a single qubit is determined by the difference between 
the resonant frequencies of the spin to be rotated and the others (its chemical shift). A very short 
pulse having a wide excitation profile in the frequency domain affects several spins at the same 
time if their chemical shifts are small. On the other hand, the duration of the Ising gate (two-qubit 
gate) depends directly on the strength of the J-coupling constants Jj k . In our setting the chemical 
shifts values impose pulse durations of the order of 1 ms, and the J-couplings impose interaction 
periods of the order of 10 ms, restricting the pulse sequences to a maximum of approximately 
1000 single-qubit rotations and 100 two-qubit (Ising) gates. 

Designing a pulse sequence to implement exactly the desired unitary transformation would re- 
quire very long refocusing schemes to cancel out all the unwanted naturally occurring J-couplings. 
Then, the overall duration of the pulse sequence increases and decoherence effects could destroy 
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our signal. Therefore, we need to find the best trade-off between the ideal 12911 accuracy of the 
pulse sequence and its duration, and neglect small couplings. For this purpose, we used an efficient 
pulse sequence compiler to perform the phase tracking calculations and to numerically optimize 
the delays between pulses, in order to minimize the error that we introduce into the quantum 
computation by neglecting small couplings. 

We now describe the parts of the pulse sequence corresponding to the three basic steps of the 
quantum simulation. 

a. Pseudo-pure state preparation: Initially, the state of the nuclei of the trans-crotonic acid 
molecules in solution is given by the thermal distribution (Sec. El). Using the methods described 
in Ref. H we have prepared the labeled pseudo-pure state (lpp) p lpp = i c 4i C3 l C2 a^l M l H2 l Hl , 
where 1 = I — a z (i.e., 1 = |1)(1|) and = I + a z (i.e., = |0)(0|). As we will see, the state p ipp , 
having the spin of Ci in the a z state, is a good initial state for our purposes. 

b. Initialization: As mentioned in Sec. |IVJ we need only 3 qubits to simulate the Fano- 
Anderson model. These qubits must be well coupled to each other to decrease the duration of 
the corresponding Ising gates we apply to them. We have chosen the spin- 1/2 nucleus Ci to 
represent qubit-1 (i.e., the impurity) and the spin- 1/2 nucleus M to represent qubit-2 (i.e., the k 
mode). On the other hand, we have chosen the spin- 1/2 nucleus C 2 to be the ancilla qubit a, to 
take advantage of its strong coupling with the spin-1/2 nucleus Ci (qubit-1). Since the rest of the 
spins (C 4 , C 3 , H 2 , Hx) in the molecule remain in the state 1 or during the whole duration of the 
experiment, we need to consider only the spins C 2 <8> Ci <8> M with the above identification. 

The initial state |+) a <g> |li0 2 ) (Sec. HvJ) can be written as p\ nit = \[{P + a^^O 2 ] in terms 
of Pauli operators. The ancilla qubit is only a control qubit and its state (i.e., its reduced density 
matrix) becomes correlated with the rest of the qubits. Since the identity part is not observable, we 
considered p inA = cr^.l 1 2 instead of p' inh as the initial state. Its preparation was done by applying 
a sequence of elementary gates to p\ pp = Vail 2 , as shown in Fig.fT2l 

c. Evolution pulse sequence: As shown in Fig. [TTJ the pulse sequence used for obtain- 
ing S(t) (Eq. [T3t requires Ising gates with a coupling strength depending on t. The refocusing 
schemes are then optimized differently and the results for different values of t cannot be directly 
comparable. To avoid this problem we have replaced the two Ising gates by an equivalent sequence 
of elementary gates, where the dependence on the simulation parameter t is transferred into the 
angle of a single-qubit rotation along the z-axis (Fig. fT3l . This virtual rotation is implemented 
through a phase tracking, as mentioned in Sec. [n] Thus, the only difference between the pulse se- 
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quence used to measure S(t) for different simulation times U is a phase calculation that introduces 
no extra optimization or experimental error. 

d. Measurement: The result of the algorithm is encoded in the polarization of the ancilla 
qubit (2er+) = (a a x ) + i(cTy) (Sec. Hill), which is directly proportional to the polarization of C 2 over 
the sample. This component precesses at the C 2 Larmor frequency u C2 - To measure it, we have to 
perform a Fourier transformation on the measured FID and integrate only the peak located at u C2 ■ 
Nevertheless, the absolute value of this signal is irrelevant since it depends on many experimental 
parameters such as the solution concentration, the probe sensitivity, and the gain of the amplifier. 
The relevant quantity is its intensity relative to a reference signal given by the observation of the 
initial state p init . To get a good signal-to-noise ratio, each experiment (or scan) was done several 
times and the corresponding experimental data were added. 

Moreover, to average over small magnetic fluctuations occurring within the duration of the 
whole experiment we interlaced scans of the reference experiment (i.e., the measurement of the 
reference signal) with scans of the actual complete pulse sequence. To increase the spatial homo- 
geneity of the field over the sample we also have inserted several automated shimming periods 
consisting of fine tuning of small additional coils located around the sample. 

B. Results 

Correlation function: In the first experiment we measured the correlation function G(t) (Eq. 
|9]) for two different sets of parameters in the Hamiltonian of Eq. |7] e fco = —2, e = —8, V = 4, 
varying t from 0.1 s to 1.5 s using increments of At = 0.1 s, and e ko = —2, e = 0, V = 4, varying 
t from 0.1 s to 3.1 s with At = 0.1 s. The duration of the optimized pulse sequences from the 
beginning of the initialization step to the beginning of the data acquisition, was 97 ms. In Fig. 
we show the analytical form of G(t) as well as the simulated and experimental data points. 
The simulated data points were obtained by a numerical simulation of the Hamiltonian dynamics 
of the full seven-qubit register under the optimized pulse sequence. This simulation is of course 
inefficient but still tractable on a conventional desktop computer. 

Hamiltonian spectrum: In the second experiment we measured the function S(t) of Eq. fl~5lto 
determine the eigenvalues of Eq. Q for Sk = —2, e = —8, and V = 0.5. The pulse sequence 
applied is the one corresponding to the quantum circuit shown in Fig. with the corresponding 
refocusing pulses. Its duration was about 65 ms. We have repeated this experiment for 128 differ- 
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ent values of the parameter t (Eq. \T5l . from t = 0.1 s to 12.8 s, using increments of At = 0.1 s. 

In Fig. [15] we show the analytical, numerically simulated, and experimental results for the 
evaluation of S(t). As mentioned in Sec. [Oil a DFT needs to be performed in order to extract the 
corresponding eigenvalues. In Fig. [T^] we show the DFT of the experimental data (see Appendix 
lAl . which reveals the expected peaks at the frequency of the two eigenvalues of Eq. |7]in the 
one-particle sector, for the above parameters. 

Discussion: At the experimental points, the error bars depend directly on the signal-to-noise 
ratio of our experimental data, as it is obtained after a fit to the experimental measured FID. They 
can then be reduced simply by running more scans for each experiment. All presented results have 
been obtained after 8 scans. 

Two different classes of errors affect the accuracy of the experimental results. The first, purely 
experimental, type of error is due to the finite accuracy of the spectrometer, and the intrinsic 
decoherence of the physical system we are working with. The second type of error is due to the 
incomplete refocusing induced by the numerical optimization scheme we used to optimize the 
pulse sequence. The numerical simulation of the optimized pulse sequence includes the errors of 
the second class but does not take into account the purely experimental ones. Thus, in our case, the 
good agreement between experimental results and simulations suggests that the main contribution 
to errors comes from the incomplete refocusing in the optimization procedure. Increasing the 
number of refocusing pulses might have led to more accurate results but would have increased the 
overall duration of the pulse sequences. The good agreement between experiment and simulation 
is consistent with the fact that the current duration of the pulse sequences are much smaller than 
the relevant relaxation time of the system (T 2 *). 

VI. CONCLUSIONS 

We have successfully simulated a quantum many-fermion system using a liquid-state NMR 
based QIP. The algebraic mapping of the operators describing any anyonic system onto the Pauli 
operators describing our QIP, combined with indirect measurement techniques, allow us to design 
efficient algorithms to simulate arbitrary evolutions of many -body anyonic systems. 

In this work the system studied was the fermionic Fano-Anderson model, which can be mapped 
onto a two-qubit system by use of the standard Jordan-Wigner transformation. Relevant dynamical 
correlation functions of the form G(t) = {4>\T^ A{T B j\<p) can be obtained by executing quantum 
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algorithms based on indirect quantum measurements, i.e., using an additional ancilla qubit. Then, 
the algorithm needed to simulate this particular system requires three qubits. We were able to 
design and run pulse sequences to implement those algorithms on an NMR QIP based on the 
trans-crotonic acid molecule (a seven-qubit quantum register). The results obtained agree with 
the theoretical ones within efficiently controlled errors. To keep a constant error level, each pulse 
sequence has been transformed such that the time parameters enter as a phase dependence. To 
shorten the duration of the pulse sequence and decrease the effect of decoherence we used only an 
approximate refocusing scheme. We numerically optimized those pulse sequences to minimize the 
error they introduce in the quantum simulation. These techniques allowed us to get very accurate 
results with efficiently controlled errors, since the overall duration of the pulse sequence was much 
smaller than the decoherence time of the system. 

Although the addition of particle-particle (e.g., density-density or exchange) interactions in 
the Fano-Anderson Hamiltonian makes it, in general, non-integrable, the quantum simulation of 
G(t) remains efficient, i.e., with polynomial complexity. We can therefore conclude that this work 
constitutes an experimental proof of principle for efficient methods to simulate quantum many- 
body systems with quantum computers. 

We thank J. Gubernatis for useful discussions on this subject. Contributions to this work by 
NIST, an agency of the US government, are not subject to copyright laws. 

APPENDIX A: DISCRETE FOURIER TRANSFORM AND PROPAGATION OF ERRORS 

Theoretically, the function S(t) of Eq. [15] is a linear combination of two complex functions 
having different frequencies: S(t) = |7i| 2 e~ iAl * + |72| 2 e~* A2 *, where A; are the eigenvalues of the 
one-particle eigenstates, defined as |lPj), in the Fano-Anderson model with n — 1 site and the 
impurity (see Sec. HV|), and A, = |(0|1P;>| 2 (Sec. M, with |0) = | X i Ts> Hi- However, the liquid 
NMR setting used to measure S(t) experimentally adds a set of errors that cannot be controlled, 
and the function S(t) shown in Fig. [15] is no longer a contribution of two different frequencies 
only. 

As mentioned in Sec. IV BL S(i) was obtained experimentally for a discrete set of values tj = 
jAt, with j = [!,■■■ ,M — 128] and At = 0.1 s. Its DFT is given by 
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where S(tj) is the experimental value of S(t) at time tj, and r/i = (with / = [1, ■ ■ • , M]) 
are the discrete set of frequencies that contribute to S(t) [I30fl . Notice that since we are evaluating 
the spectrum of a physical (Hermitian) Hamiltonian, the imaginary part of S(rji) is zero 13 ill . In 
Fig- El we show S(rji) obtained from the experimental points S(tj) of Fig. [l5j Its error bars (i.e., 
the size of the line in the figure) were calculated by considering the experimental error bars of 
S(tj) in the following way: First, we rewrite Eq. lAll as 



S{vi 



M 
3=1 



(A2) 



with Qij = M 1 [Re(S(tj)) cos(r]itj) — \m(S(tj)) sm(r]itj)] (real). Then, the approximate standard 
deviation ES t of S(rji) depends on the errors EQij of Qij as (considering a normal distribution 

1220) 

M 

[EStf » ^[EQ^] 2 . (A3) 
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On the other hand, EQu is calculated as 113211 



dQij 



dRe(S(tj)) 



Er 2 + 



dQij 



d\m(S(tj)) 
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where E R and E| are the standard deviations of the real and imaginary parts of S(tj) (see Fig.fT3t. 



respectively. Because of experimental reasons (Sec. IV Al) these errors are almost constant, having 
E R ~ E| ~ Es independently of tj (see Fig. IT5t . where E$ is taken as the largest standard deviation. 
Combining Eqs. IA3l and lA4l we obtain 



ES, 



M~ 



M 

! E s 2 J][|cos(77^)| 2 + |sm(77^; 

3=1 



1/2 



M 



(A5) 



In our experiment, M = 128 and E s ~ 0.04, obtaining E5/ ~ 0.0035, which determines the 
(constant) error bars (i.e., the size of the dots representing data points) shown in Fig.lTol 

The standard deviation Er]i in frequency domain is due to the resolution of the sampling time 
At. This resolution is related to the error coming from the implementation of the ^-rotations in the 
refocusing procedure (Fig. El). A bound for this error is given by the resolution of the spectrum; 
that is, 
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FIG. 1: Circuit representation of the elementary gates. The top picture indicates a single-qubit rotation 
while the bottom one indicates the two-qubit Ising gate. Any quantum algorithm can be represented by a 
circuit composed of these elementary gates (see for example Fig.0 
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II) 

FIG. 2: Bloch's sphere representation of a single nuclear spin-1/2 precessing around the quantization axis 
determined by the external magnetic field B. The precession frequency is given by Uj = fJ,jB, with /ij the 
magnetic moment of the j-th nucleus. Due to the chemical environment, each nucleus precesses at its own 
Larmor frequency Uj. 
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FIG. 3: Circuit representation for the refocusing scheme to control J-couplings. The Ising-like coupling Jjk 
between spins can be controlled by performing flips on one of the spins at times t\ = At\ and t% = t±+ At2, 
respectively. The effective coupling is a = a± — 02 = Jjk(At± — At2), and vanishes when At± = At2- 
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(2al) = (4>\U(t) 



U(t) 



FIG. 4: Quantum network for the evaluation of the expectation value of a unitary operator U (t). The filled 
circle denotes a controlled operation (i.e., U' 1 ^ of Sec. HHt . such that U{t) is applied to the system only if 
the ancilla qubit is in the state |l) a . 
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FIG. 5: Quantum network for the evaluation of the spectrum of an observable Q. 
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FIG. 6: Quantum network for the evaluation of the correlation function G(t) = {<fi\T^ A(I 'B j\4>) . The filled 
(empty) circle denotes an operation controlled in the state |l) a (|0) a ) of the ancilla qubit. 
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FIG. 7: Mapping of the fermionic product state c[c 2 c\\vac) , with |vac) the no-fermion or vacuum state, into 
the spin- 1/2 and the standard quantum computation languages, using the Jordan-Wigner transformation. A 
filled circle denotes a site occupied by a spinless fermion, which maps into the state |f) in the spin 1/2 
algebra. 
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FIG. 8: Fermionic Fano- Anderson model. Fermions can hop between nearest-neighbor sites (exterior cir- 
cles) and between a site and the impurity (centered circle), with hopping matrix elements r and V/y/n, 
respectively. The energy of a fermion in the impurity is e. 
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FIG. 9: Quantum circuit for the evaluation of (Eq. |9j in terms of elementary gates directly imple- 
mentable with liquid-state NMR methods. 
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FIG. 10: Quantum circuit for the evaluation of S(t) (Eg. [T3t. The parameters Ai and A2 are defined in Sec. 
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. The decomposition of the operator U in NMR gates can be found in Fig. |9] 
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FIG. 1 1 : The trans-crotonic acid molecule is a seven-qubit register: The methyl group is used as a single 
spin 

1/2 fl 

and four 13 C. The table shows in hertz the values of the chemical shifts (on the main diagonal) 
and the J-couplings (off-diagonal) between every pair of nuclei (qubits). 
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FIG. 12: Initialization pulse sequence used to transform the initial labeled pseudo-pure state p\ pp = 
1 2 a^l M into the state p\ n - lt = <r£ 2 l Cl M . The sequence transfers the polarization from Ci to C2 and 
flips the spin of the methyl group M. We have chosen the spin- 1/2 nuclei C2, Ci, and M to represent the 
ancilla, qubit-1 (i.e., the impurity), and qubit-2 (i.e., the /co-mode), respectively. 



30 




FIG. 13: Modification of a two-qubit gate with a coupling strength depending on a parameter t. The variable 
interaction period is translated into fixed interaction periods and a single-qubit rotation with variable angle 
about the z-axis. Using this trick, the duration of the physical pulse sequence does not depend on the 
parameter t representing the time of the simulation. 
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FIG. 14: Real and imaginary parts of the correlation function G(t) of Eq.|9] The top panels show the results 
when the parameters in Eq. 0are £fc = —2, e = —8, V = 4. The corresponding parameters Ax, A2, are 
in the quantum network, Fig. ^] are used to measure G(t) and can be determined using Eqs. [T3land[T4l 
The bottom panels show the results for ek = —2, e = 0, V = 4. The (black) solid line is the analytic 
solution, the red circles are obtained by the numerical simulation (including the refocusing pulses), and the 
blue circles with the error bars are experimental data. 
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FIG. 15: Real and imaginary parts of S(t), for ek = —2, e = —8, and V = 0.5 in Eq. The (black) solid 
line corresponds to the analytic solution. The red circles correspond to the numerical simulation (using 
refocusing pulses) and the blue circles with the error bars are experimental data. S(t) has been measured 
using the network of Fig. [TO] with a = (e + Sk )/2. 
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FIG. 16: Discrete Fourier transform of the real part of the experimental data of Fig. \^\ The position of 
the two peaks corresponds to the two eigenvalues of the Hamiltonian of Eq. 0for Ek = — 2, e = —8, and 
V = 0.5. Numbers in parentheses denote the exact solution. The size of the dots representing experimental 
points is the error bar (see Appendix©. An upper bound to the error in the frequency domain is 0.5, 
which was determined by the resolution of the spectrum. 
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